﻿using MathNet.Symbolics;
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;

using MathNet.Numerics.LinearAlgebra;
using MathNet.Numerics.LinearAlgebra.Double;
using MathNet.Numerics.Statistics;
using MathNet.Symbolics;
using Expr = MathNet.Symbolics.SymbolicExpression;
using MathNet.Numerics.Integration;
using System.Linq.Expressions;
using static System.Runtime.InteropServices.JavaScript.JSType;


namespace NumericalAnalysis.DataInterpolation
{


    public class HermiteInterpolation : IInterpolation
    {
        static Expr polynomial = 0;   //多项式表达式
        List<double> polyCoefficients = new List<double>(); //各项系数 幂次高至低
        List<int> coefficienOrders = new List<int>(); //各项幂次

        static Vector<double> theVec = Vector<double>.Build.Dense(2);
        List<Vector<double>> vectors = new List<Vector<double>>();

        List<double> xs = new List<double>();
        List<double> ys = new List<double>();
        List<double> yds = new List<double>();

        int count = 0;

        static string var_t { get; set; } = "t";

        //static string IInterpolation.var_t { get; set; } = "t";

        /// <summary>
        /// 初始化
        /// </summary>
        /// <param name="_vectors"></param>
        public HermiteInterpolation(IEnumerable<Vector<double>> _vectors, IEnumerable<double> _yds)
        {
            count = _vectors.Count();
            vectors = _vectors.ToList();

            xs = _vectors.Select(v => v[0]).ToList();
            ys = _vectors.Select(v => v[1]).ToList();

            if (_yds.Count() == count)
            {
                yds = _yds.ToList();
            }
            else
            {
                for (int i = 0; i < count - 1; i++)
                {
                    double diff = (ys[i + 1] - ys[i]) / (xs[i + 1] - xs[i]);
                    yds.Add(diff);
                }
                yds.Add(yds.Last());
            }
        }

        public HermiteInterpolation(IEnumerable<Tuple<double, double>> _tuples, IEnumerable<double> _yds)
        {
            count = _tuples.Count();
            vectors = _tuples.Select(v =>
            {
                var vec = Vector<double>.Build.Dense(2);
                vec[0] = v.Item1;
                vec[1] = v.Item2;
                return vec;
            }).ToList();

            xs = _tuples.Select(v => v.Item1).ToList();
            ys = _tuples.Select(v => v.Item2).ToList();

            if (_yds.Count() == count)
            {
                yds = _yds.ToList();
            }
            else
            {
                for (int i = 0; i < count - 1; i++)
                {
                    double diff = (ys[i + 1] - ys[i]) / (xs[i + 1] - xs[i]);
                    yds.Add(diff);
                }
                yds.Add(yds.Last());
            }
        }


        /// <summary>
        /// 常规公式解法
        /// </summary>
        /// <returns></returns>
        public Expr FitInterpolation()
        {
            //polynomial = 0;
            if (polynomial.ToString().Contains(var_t)) return polynomial;

            var t = Expr.Variable(var_t);

            for (int i = 0; i < count; i++)
            {

                Expr h_i = 1;
                double a_i = 0;
                for (int j = 0; j < i; j++)
                {
                    h_i *= ((t - xs[j]) / (xs[i] - xs[j])).Pow(2)  ;
                    a_i += 1 / (xs[i] - xs[j]);
                }
                for (int j = i + 1; j < count; j++)
                {
                    h_i *= ((t - xs[j]) / (xs[i] - xs[j])).Pow(2);
                    a_i += 1 / (xs[i] - xs[j]);
                }

                double x_i = xs[i];
                double y_i = ys[i];
                double yd_i = yds[i];
                //yd_i= polynomial.DifferentiateAt(t,x_i);

                Expr basisFunc = h_i * ((x_i - t) * (2 * a_i * y_i - yd_i) + y_i);

                polynomial += basisFunc.Expand();
            }

            //显示多项式
            Console.WriteLine("--------------------------------");
            Console.WriteLine("多项式 polynomial:");
            Console.WriteLine(polynomial);

            //Console.WriteLine("--------------------------------");
            //Console.WriteLine("多项式简化 polynomial-ExponentialSimplify:");
            //Console.WriteLine(polynomial.ExponentialSimplify());

            //Console.WriteLine("--------------------------------");
            //Console.WriteLine("多项式展开式 polynomial-Expand:");
            //Console.WriteLine(polynomial.Expand());

            polynomial = polynomial.Expand();

            return polynomial;
        }


        /// <summary>
        /// 从给出的向量计算
        /// </summary>
        /// <param name="sourceVecs"></param>
        /// <returns></returns>
        public IList<Vector<double>> CalculateBy(IList<Vector<double>> sourceVecs)
        {
            for (int i = 0; i < sourceVecs.Count(); i++)
            {
                double xValue = sourceVecs[i][0];
                var symbols = new Dictionary<string, FloatingPoint> { { "t", xValue } };
                double yValue = polynomial.Evaluate(symbols).RealValue;
                //double   yValue = CalculateBy(xValue);
                sourceVecs[i][1] = yValue;
            }
            return sourceVecs;
        }

        /// <summary>
        /// 从给出的数列计算
        /// </summary>
        /// <param name="xValues"></param>
        /// <returns></returns>
        public IList<double> CalculateBy(IList<double> xValues)
        {
            List<double> yValues = new List<double>();
            for (int i = 0; i < xValues.Count(); i++)
            {
                var symbols = new Dictionary<string, FloatingPoint> { { "t", xValues[i] } };
                double yValue = polynomial.Evaluate(symbols).RealValue;

                yValues.Add(yValue);
            }
            return yValues;
        }

        /// <summary>
        /// 从给出的数计算
        /// </summary>
        /// <param name="xValue"></param>
        /// <returns></returns>
        private double CalculateBy(double xValue)
        {
            //复杂多项式避免使用Compile方法计算  将严重拖慢性能
            //Func<double, double> func = polynomial.Compile("t"); 
            //double yValue = func(xValue);

            var symbols = new Dictionary<string, FloatingPoint> { { "t", xValue } };
            double yValue = polynomial.Evaluate(symbols).RealValue;

            return yValue;
        }


    }

    public class HermiteInterpolationTest 
    {
        string var_t = "t";

        /// <summary>
        /// sin函数
        /// </summary>
        /// <returns></returns>
        public (List<Tuple<double, double>>, List<Tuple<double, double>>) ExampleSin()
        {
            List<double> interpXs = NumericalService.LinSpace(0, Math.PI * 2, Math.PI * 2 * 0.1);
            List<double> checkXs = NumericalService.LinSpace(0, Math.PI * 2, Math.PI * 2 * 0.02);   //构造sin函数插值点 

            var t = Expr.Variable(var_t);
            Expr sourceExpr = t.Sin();  //原函数
             
            //原函数的离散点
            var interpPnts = sourceExpr.GeneratePointBy(var_t, interpXs).ToList();

            ////进行插值  
            Expr diffExpr = sourceExpr.DiffBy(var_t);
            var yds = diffExpr.GeneratePointBy(var_t, interpXs).Select(d=>d.Item2);
            //yds = new List<double>(); ;
            var interp = new HermiteInterpolation(interpPnts , yds);
            var expression = interp.FitInterpolation();  //插值函数
            var checkdPnts = expression.GeneratePointBy(var_t, checkXs); // interpTest.GenerateCheckPoints(interp, array2);

            return (interpPnts.ToList(), checkdPnts.ToList());
        } 

        public (List<Tuple<double, double>>, List<Tuple<double, double>>) ExampleRunge()
        { 
            List<double> interpXs = NumericalService.LinSpace(-5, 5, 0.5);
            List<double> checkXs = NumericalService.LinSpace(-4, 4, 0.2);
           
            var t = Expr.Variable(var_t); 
            Expr sourceExpr = "1 / (t * t + 1)";  //原函数

            //原函数的离散点
            var interpPnts = sourceExpr.GeneratePointBy(var_t, interpXs).ToList();

            ////进行插值  
            Expr diffExpr = sourceExpr.DiffBy(var_t);
            var yds = diffExpr.GeneratePointBy(var_t, interpXs).Select(d => d.Item2);
            var interp = new HermiteInterpolation(interpPnts, yds);
            var expression = interp.FitInterpolation();  //插值函数
            var checkdPnts = expression.GeneratePointBy(var_t, checkXs); // interpTest.GenerateCheckPoints(interp, array2);

            return (interpPnts.ToList(), checkdPnts.ToList());

        }

        public (List<Tuple<double, double>>, List<Tuple<double, double>>) ExampleHomeWork()
        {
          
            List<Tuple<double, double>> interpPnts = new List<Tuple<double, double>>();
            List<double> values = new List<double>()
                { 12, 9, 9, 10, 18, 24, 28, 27, 25, 20, 18, 15, 13 };
            for (int i = 0; i < values.Count; i++)
            {
                double xk = i * 2;
                double yk = values[i];
                Tuple<double, double> tuple = new Tuple<double, double>(xk, yk);
                interpPnts.Add(tuple);
            }
              
            List<double> checkXs = NumericalService.LinSpace(0 + 0.1, 24 - 0.1, 0.1);
             
            ////进行插值   
            var interp = new HermiteInterpolation(interpPnts, new  List<double>());
            var expression = interp.FitInterpolation();  //插值函数
            var checkdPnts = expression.GeneratePointBy(var_t, checkXs); // interpTest.GenerateCheckPoints(interp, array2);

            return (interpPnts.ToList(), checkdPnts.ToList());

        }

    }
       
}
